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Abstract 

Subgrid analysis of a transitional temporal mixing layer 
with evaporating droplets has been performed using a di- 
rect numerical simulation (DNS) database. The DNS is 
for a Reynolds number (based on initial vorticity thick- 
ness) of 600, with droplet mass loading of 0.2. The gas 
phase is computed using a Eulerian formulation, with La- 
gr&ngian droplet tracking. Since Large Eddy Simulation 
(LES) of this flow requires the computation of unfiltered 
gas-phase variables at droplet locations from filtered gas- 
phase variables at the grid points, it is proposed to model 
these by assuming the gas-phase variables to be given by 
the filtered variables plus a correction based on the filtered 
standard deviation, which can be computed from the sub- 
grid scale (SGS) standard deviation. This model predicts 
the unfiltered variables at droplet locations better than 
simply interpolating the filtered variables. Three methods 
are investigated for modeling the SGS standard deviation: 
Smagorinsky, gradient and scale-similarity. When prop- 
erly calibrated, the gradient and scale-similarity methods 
give results in excellent agreement with the DNS. 

Introduction 

Droplet-laden turbulent flows occur in many problems of 
practical interest, such as spray combustion and atomiza- 
tion. The interaction of particles and turbulence is an inte- 
gral feature of such flows, and hence the topic of much re- 
search [2] [3] [5] [12] . Large Eddy Simulation (LES), in which 
the flow field is spatially filtered, is emerging as a powerful 
tool in modelling unsteady turbulent flows. It is expected 
to be more generally applicable than Reynolds- Averaged 
Navier Stokes (RANS), since the large scale structures 
are computed, and the more universal small scale struc- 
tures are modelled. LES is also less computationally in- 
tensive than direct numerical simulation ( DNS) in which 
ail length scales are resolved, and has the additional ad- 
vantage of being able to accommodate considerably larger 
Reynolds numbers (Re). Whereas much research has been 
devoted to LES modeling for single phase incompressible 
flows, only moderate attention has been given to compress- 
ible shear flows [1] [15], with focus now turning two-phase 
flows [4] [11] [16]. In addition to modeling subgrid scale 
(SGS) terms for the gas phase, an LES of a droplet-laden 
flow would require modeling the unfiltered gas phase vari- 
ables at the droplet locations. In simplistic models, the 


filtered variables are substituted for the unfiltered vari- 
ables; however, this assumption may be substantially in- 
accurate for droplets with small Stokes numbers. With an 
increasing body of DNS computations [5] [6] [7] [8] [13] [14], it 
is now possible to assess subgrid scale quantities at mod- 
erate Reynolds numbers, with good prospects for devising 
subgrid scale models. 

Recently, Miller and Bellan [9] [10] have performed 
DNS of droplet laden mixing layers. They use ‘DNS’ to 
refer to computations in which all length scales of the 
gas-phase are resolved but the effect of the gas on each 
droplet is modeled using Stokes drag, and the effect of 
the droplets on the gas are modeled as source terms in 
the gas-phase equations. The present paper addresses the 
use of the DNS database of Miller and Bellan [9] to eval- 
uate subgrid scale closures. Specifically, we examine the 
largest Reynolds number (based on initial vorticity thick- 
ness, 6 Wt o, and initial velocity difference, AC/ o) of 600, 
with mass loading of 0.2 (3x 10 6 drops) on a 300 x 332 x 180 
(0.25m x 0.22m x 0. 12m) grid, which they denote as case 
TP600. 


Governing Equations 


The governing equations for the gas phase density (/?), ve- 
locity (ui), total energy ( E ) and vapor mass fraction (Yy) 
are given by: 
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where subscript V denotes the vapor, subscript C denotes 
the carrier g;us, and the mass fraction of the carrier gas is 

Y c = l - Y v . 

The Lagrangian particle equations for the position 
(X t ), velocity (u,) temperature (T^) and mass (m^) are: 
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Computation of the drag force F», the heat flux Q and the 
evaporation rate require knowledge of the gas phase 
variables (u iy T t Yv,P) at the droplet locations, and in- 
volve validated relations [10]. 

The source terms are 
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where 51 indicates appropriately weighted summations 
over droplets within a discretization volume associated 
with each grid point. 


Filtered Governing Equations 
The filtering operation is defined as: 

0(“?) = [ - ~y )dy* 
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where G is the filter function, with the property that 1 = 1, 
and V is the filtering volume. We use a cubic top-hat 
filter, in which V is a cube of sides A, and G is simply 
a volume-average. For compressible flow, we use Favre 
filtering, defined as ^ = p4>/]$, to simplify the notation. 
After filtering, the gas phase equations become: 
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and it has been assumed that 
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The terms that need to be modeled are the subgrid stress 
Tij , the subgrid heat flux 9j and the subgrid species flux 


The droplet equations remain as before, except that 
now the gas-phase variables (u» , T y Yy , P) at the droplet lo- 
cations are no longer immediately available, and^wiU need 
to be derived from the filtered variables (tXi,T, Yy, P). 
Thus the (unfiltered) gas-phase variables also need to be 
modeled. 


Model for Instantaneous Variables 

The droplet model requires modeling gas phase variables 
at droplet locations. To guide the modeling, we will first 
consider the known DNS generic variable <j> and its filtered 
form where the bar denotes Favre filtering for u, and 
Yy and regular filtering for T and P. The definition of the 
standard deviation is 

<T = \J = yj (<j> - <p) 2 (35) 

Thus the relation between (p and <t> is 

0 = + /<T (36) 

where from the definition of <r, / = ±1. The goal of the 
modeling is to compute, from the filtered flowfield, the 
form of fa that provides a better approximation to <f> than 
does fa = 0. In this formulation, fa can be viewed as a 
correction to ^with sign / and magnitude a. 

It is tempting to assume that / randomly takes on 
values of -1 and 1. However, if the filtering operation is 
viewed as a volume average, a relation between <p and <p 
can be derived as follows. Consider the third-order Taylor 



expansion of </> in the filtering volume V of size A with 
centroid at r 0 - (.c lt , , x 2 „ , xrj„ ) integrated over the volume: 
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where terms of 0(A 3 ) vanish due to the symmetry of the 
filtering volume. 

Thus / will generally be —sign (V 2 tf>). From the 
available filtered quantities, we can compute V 2 <£ rather 
than V 2 0; so to model / we assume that V 2 <£ and V 2 <£ 
have the same sign. To model a*, we note that for the gas 
phase we will be modeling terms of the form <fxf> — <j> <f>, 
which are the subgrid scale fluctuations. Defining <tsgs 83 
the SGS standard deviation, 
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If we assume that R{<t>, <f>) — 1, then 



Defining rr = , and using 7 as a model for <r, we 

arrive at a model for <fi of the form 


<p = <f> — sign(V 2 <t>y& (46) 


i.e. / = -aiyn(V 2 0), rr = 7f. To assess this model, we will 
use the “exact” 7f obtained by filtering the DNS flowfield, 
and then turn our attention to modeling it based on the 
filtered flowfield. The proposed model will be computed 
at the grid points, and then interpolated to the droplet 
locations. 

Figure 1 shows the probability density function(PDF) 
for {<(> — <i>) ft? for case TP600 at the end of the simulation. 
These results were obtained using a cubic top-hat (box) 
filter with filter width A = 4Ax. Notably, for all quantities 
there are peaks near 1, and for the velocity components 
there are also peaks near -1, similar to the PDF of (<f> — 
0)/<t, which takes on values of ±1. 




Figure 1: PDF of (4> — 4>)/rr 

Figures 2 and 3 show the comparison between inter- 
polating the DNS flowfield to the droplet locations (the 
“exact” quantities) and interpolating the models to the 
droplet locations. Results are presented in terras of av- 





nnig<*s over droplets within a given y-intervnl; this aver- 
aging is denoted by <<>>. First, using rr os given by 
Equation 35, we compare the models for /. It is seen 
that / = 0 leads to significant discrepancy between the 
exact and model interpolated variable. Any model with 
/ with mean 0 will not perform any better as the devise 
tions toward the exact value will be just as likely as the 
deviations away from the exact values. This is seen in the 
case where / is randomly taken to be -1 or +1, which per- 
forms slightly worse than / = 0 despite having the exact 
rr. Using / = -si0n(V 2 0) gives significant improvement 
in view of the two assumptions made i.e. / = -sign(V 2 <t>) 
and sign(V 2 <t>) = si^n(V 2 0). For the data shown in Fig- 
ures 2 and 3, this is true (86%, 86%, 86%, 85%, 86% and 
93%) of the time for ui,ix 2 ,u 3 ,7\ W, and P respectively. 
When rr is replaced by <r (denoted [cr] in the figures), with 
/ = -siyn(V 2 0), there is considerable improvement over 
/ = 0, for all quantities except T, where all the models give 
similar results. For T, Yy and P, / = — sign( V 2 <£) can be 

replaced by / = 1 if the signed <7 = y/^4 — \J 4>4> is used 

instead of the positive 7? = \f<P4> — \f^4> of Equation 

45. An alternative expression for / uses a £cale*similarity 
idea, i.e. / = sign{<}> — <j>) = sxgn(<t>_~ <p). This gives 
similar results to using / = — sign(V 2 (p) for all six vari- 
ables. An analysis of the correlation R(<p, (p) of Equation 
44 showed that R = 1 for T and P, 0.97 < R < 1 for ui, u 2 
and u 3 and 0.7 < R < 1 for Yy. For t*i,u 2 and u 3 , the 
greatest deviation from 1 was in the middle of the mixing 
layer. For Yy , the greatest deviation was at the droplet- 
laden/droplet-free interface. However, even in this region, 
R = 1 led to the best prediction of Yy at droplet locations 
compared to smaller values. Most importantly, R = 1 
is the only value that provides the correct <f ^ 0 in the 
laminar freestreara. Finally, we note that the error for the 
approximation of Equation 46 is determined by that of the 
velocity components, which have the largest error of about 
1.5%. 


Models for Subgrid Cross-term s 

For LES in the gas-phase, models are required for the 
subgrid stresses Tij = — u%u heat fluxes 9j = 

Tuj - Tuj and species fluxes g 3 = YyUj - Yyiij. For the 
droplet part, models are required for the subgri d va riances 
(T 2 S cs : - UlUi, — U 2 U 2 , U3U3 ~ U 3 U 3 , TT - TT, 

YyYy — YyYy and PP-P P. Since these two sets of terms 
are of the same form, it seems reasonable and consistent 
to use the same type of model for both. We consider three 
models, Smagorinsky, Gradient, and Scale-Similarity [15]. 
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Figure 2: Error in unfiltered variable model (p = <P + f<* 
interpolated to droplet locations 







Figure 3: Error in unfiltered variable model 4> = <t> + fa 
interpolated to droplet locations 


1, = ~Cr tfy/SuSu^- (49) 

with model constant Cn and filter width A, and where the 
rate-of-strain tensor is defined as 
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Although this is forms the basis for most SGS modeling, 
it cannot be easily extended to compute the subgrid vari- 
ances for T, P and Yy . 


Gradient SGS Model 
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This model is easily extended to compute the subgrid vari- 
ances for any quantity <f> as 

°sgs W = 4><t> — <$><$> — CgA 2 — — (54) 

Theoretically, CgA 2 is the moment of the filtering volume, 
I c A 2 of Equation 40, as can be seen by integrating the 
square of the Taylor expansion for 0, Equation 38, over 
the filtering volume, but using filtered quantities in the 
derivatives. Thus for a cubic top-hat filter Cg = ^ * This 
model has the advantage that the derivatives are already 
available from the computation of the resolved part. 

Scale-Similarity SGS Model 

This model involves refiltering the flow-field at a test filter 
A > A _ 

Tij = Cs (uiUj - (55) 

6, = C s (^f- ufi) (56) 

r 1j = C s (tjv ~ %Y v ^j (57) 

Theoretically, Cs should be 1 . This model is easily ex- 
tended to compute the subgrid variances for any quantity 
<p as 

^scs (0) = 4><t> — <P<P = Cs {&<t> - <t> (58) 

Model Coefficients 

The validation of the models involves comparing the values 
predicted by the models to those obtained from the DNS 
database. Standard deviations from Equation 42 will be 
referred to as “exact”. First, the correlation between the 
“exact” and model standard deviations will be computed. 






Variable 

Slope= y/Ca 

Correlation 

'M 

0.40.17 

0.9855 

u 2 

0.4087 

0.9847 


0.4112 

0.9835 

T 

0.3886 

0.8727 

Yv "1 

0.4155 

0.9806 

P 

0.3980 

0.9902 

Combined 

0.40 



Table 1: Gradient SGS Model A = 4Az 


The correlations are computed by averaging over homoge- 
neous (ii — xz) planes 
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or over the whole domain 
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By definition R is between -1 and 1. Values near 
1 indicate strong positive correlation, values near -1 in- 
dicate strong negative correlation, whereas values near 0 
indicate poor correlation. This allows pointwise assess- 
ment of the correlations. Next, the relationship between 
the two variables (i.e. the model coefficient) needs to be 
determined. The simplest is a constant-coefficient model 
where Cr,C g , Cs are the same for all flow variables over 
the whole (spatial, temporal) domain. In that case, the 
coefficient can be determined using a least-squares fit to 
Y = bX which leads to b = [XY] /[XX]. If X is the model 
standard deviation and Y is the “exact” standard deviar 
tion , then b is the square root of the model coefficient. 
More sophisticated models would have the model coeffi- 
cients as functions of space and time. 


Results 

All three models were evaluated for Reynolds number of 
600 and mass loading 0.2 (case TP600 of Miller and Bellan 
[9]) at the end of the simulation, using cubic top-hat filters. 

Figure 4 shows the subgrid scale stresses predicted by 
the Smagorinsky model. It can be seen that there is little 
local correlation between the exact and predicted stresses. 
Thus, this model was not assessed further. 

Figures 5 and 6 show the plane- aver aged SGS stan- 
dard deviations for ui,tx 2 ,ti 3 and T,Yy,P respectively. 
Deviations from both models show good correlations with 
the exact deviations, although the models tend to under- 
predict the temperature deviations. For these figures, val- 
ues of C G = 0.4 2 and C s = 0.79 2 were used. These values 
were obtained from linear fits of the SGS standard deviar 
tions, which produced the coefficients presented in Tables 
1 and 2. 





Figure 4: Subgrid stresses a)exact b) Smagorinsky c) cor- 
relation 



















Variable 

Slop 0.= y/Cs 

Correlation 


0.7465 

0.9475 

u 2 

0.7772 

0.9498 

W:i 

0.8008 

0.9498 

T 

0.9143 

0.8958 

Yv 

0.7856 

0.9181 

P 

0.7350 

0.9476 

Combined 

0.79 



Table 2: Scale-Similarity Model A = 4Ax; A = 8Ax 
Conclusions 

An a priori subgrid analysis has been conducted for a tem- 
porally developing mixing layer with evaporating droplets. 
This analysis was performed on a DNS database for 
Reynolds number (based on vorticity thickness) of 600 and 
mass loading of 0.2. Two models for the subgrid scale 
(SGS) standard deviations, the gradient and scale simi- 
larity models, where found to give excellent results when 
the model constant was properly calibrated. A model to 
recover the unfiltered variables from the filtered variables 
was also examined. In this model, the unfiltered variables 
are taken to be filtered variables plus a correction term 
which can be computed from the SGS standard deviations. 
Predictions for the unfiltered variables at the droplet lo- 
cations were found to be improved compared to simply 
interpolating the filtered variables. Future work involves 
both testing these models on a DNS database for a higher 
mass loading and a posteriori testing of these models in an 
LES. 
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